Universal behavior in populations composed of excitable and self-oscillatory elements 
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We study the robustness of self-sustained oscillatory activity in a globally coupled ensemble of 
excitable and oscillatory units. The critical balance to achieve collective self- sustained oscillations 
is analytically established. We also report a universal scaling function for the ensemble's mean 
frequency. Our results extend the framework of the 'Aging Transition' [Phys. Rev. Lett. 93, 104101 
(2004)] including a broad class of dynamical systems potentially relevant in biology. 
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Large networks of interacting dissipative systems are 
appropriately modeled in terms of coupled nonlinear dif- 
ferential equations, which successfully reproduce a huge 
variety of the dynamical patterns found in nature. In 
many cases, the individual systems present time-periodic 
behavior and may achieve a certain degree of global syn- 
chronization despite the unavoidable differences among 
them 0, 0. This subject has attracted a great deal of 
both theoretical and experimental interest during the last 
decades Q. 

However, the robustness of the macroscopic synchro- 
nized oscillations in a mixed population of self-oscillatory 
and non-self-oscillatory elements has only been addressed 
very recently by Daido and Nakanishi 4J. Interestingly, 
they report a general scenario, called 'Aging Transition', 
characterized by a universal (i.e. independent of the os- 
cillator type) scaling function. In £| the Aging Transition 
was found in populations of oscillators in which some of 
them lose their self-oscillatory activity (by deterioration 
or 'aging') through an inverse Hopf bifurcation. 

In this Letter we present a extension of the Aging Tran- 
sition taking place in a new class of systems in which the 
oscillatory behavior is lost in a saddle-node (SN) bifurca- 
tion. Such systems are of particular relevance since their 
resulting dynamics is excitable and thus of interest in 
many areas of physics, chemistry and biology 0,|HE10- 

For instance, a remarkable example of macroscopic 
synchronization is found among the pacemaker cells in 
the sino-atrial node which initiate the heartbeat. When 
a certain ratio of cells are damaged by disease, the lack 
of an adequate synchronized state requires the implanta- 
tion of an electronic pacemaker. Nevertheless, new tech- 
niques aiming to create biological pacemakers have been 
recently proposed as a very promising alternative to elec- 
tronic ones y|. One of them consists in the creation of 
an aggregate of biological pacemaker cells in some re- 
gion(s) of the ventricle, converting excitable heart cells 
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FIG. 1: Morris-Lecar equations (J5| in the excitable (left 
panel) and in the oscillatory /pacemaker (right panel) regimes. 
Both dynamical states are linked through a saddle-node bifur- 
cation on invariant circle (SNIC) at <f>* ~ 0.076 [see Eq. I2al l, 
The nullclines V — and w — are depicted with gray lines. 

into pacemaker cells by gene transfer ■£>), in an "inverse- 
aging" transition. This is feasible since the heart's pace- 
maker and excitable cells are very similar (they are both 
identical pacemaker cells in the early embryonic heart, 
and differentiate as the development progresses). In this 
context, the robustness of the aggregate's self-oscillating 
activity in a mixed population of self-oscillatory (con- 
verted) and excitable (unaltered) elements seems to be 
particularly relevant. 

In our study, we have considered systems that cease 
their oscillatory behavior through a SN bifurcation on 
invariant circle (SNIC) @, also called SN homoclinic bi- 
furcation ^(j ( see Fig-©- This is the simplest possible 
scenario linking excitable and oscillatory dynamics, and 
there is only one attractor at each parameter value. The 
excitable regime (at one side of the bifurcation) is very 
well known in theoretical neuroscience, where it is re- 
ferred to as Class I excitability [Tl|. 

Additionally, we make the following assumptions: 

(i) The population (of size N) is divided into two 
groups, consisting of pN [1/N < p < (N — 1)/N] 
identical excitable units (Se), and of (l—p)N iden- 
tical active units (Sa)- The elements in the pop- 
ulation are ordered according to an index j, such 
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FIG. 2: (color online) (K,p) phase diagram for a mixed popu- 
lation of oscillatory and excitable elements: (a) Morris-Lecar 
@ with 4>a = 0.1 and <f>E = 0.06; (b) phase model © with 
bA — (oscillatory) and bs = \/2 (excitable). Inside the 
region of Global Oscillations the black dashed line separates 
1 : 1 (right) from other 1 : r, < r < 1 (left) frequency 
ratios, within the region GO. The red dotted curves are the 
deformation of the solid lines when heterogeneity in each sub- 
population is considered (see text): 7 = (a) 0.01 and (b) 0.6. 

that S A = j £ {1, . . . , (1 - p)N} and S E = j e 
{(l-p)N + l,...,N}. 

(ii) A linear all-to-all coupling is assumed, for the state 
variable x, of each unit: 

K N 

fe = l 

where Fj = for j G Sa and Fj = F E for j G Se- 
In general, the coupling term in Q could simply 
enter through a single variable of Xj (e.g. the mem- 
brane potential, in a cell model) . This has in many 
cases the same synchronizing effect (but sec 12]). 

First we introduce our findings using an ensemble of 
globally coupled Morris-Lecar (ML) units [Hi . The ML 
equations were obtained from the study of the electrical 
activity of the barnacle muscle fiber, and later t hey were 
popularized as a reduced model of excitability [II]]. We 
have used the adimensionalized version proposed in |14| , 
with the coupling term (proportional to K) entering 
through the voltage variable. The system is 

CVj = g L (-V L - V 3 ) + gcamooiV^Vca - Vj) (2a) 

K N 

- gKW,(V K + Vj) + 0,(0.2 - Vj) + - - Vi), 

i=l 

Wj = MVj)(woo(Vj) — Wj), (2b) 

where moo(Vj-) = 0.5{1 + tanh[(Vj — vi)/v 2 }}, Woo{Vj) — 
0.5{l+tanh[(V,-W3)M]}, andAfV,) = A {l+cosh[(I^- 
v a)/ v 4,]}- Several constants [l!j were taken from [lj] 
with the external current set equal to zero. The term 
4>j (0.2 — Vj) in Eq. I|2a|l controls the dynamics of an iso- 
lated element [2(j . In our model <pj takes only two possi- 
ble values 4>a and 4>e- For <j) A > <\>* — 0.076 the dynamics 



of an isolated cell is self-oscillatory, whereas for 4>e < 4>* 
it is excitable. Both behaviors are linked by a SNIC at 
<t>* (Fig. EE). 

The system is described in terms of two control pa- 
rameters, the ratio p of excitable units, and the coupling 
strength K. Figure^a) shows the (K,p) phase diagram 
with the dynamical regimes found in the ensemble of ML 
units. For small coupling K, each element in the popu- 
lation essentially maintains its intrinsic dynamics. Such 
state is labeled as Partial Oscillations (PO) because the 
collective state is not entirely oscillatory (the excitable 
elements are not able to oscillate, they only exhibit small 
amplitude "pulsations"). For large values of the coupling 
strength K, all the elements in the ensemble exhibit the 
same behavior, but depending on the ratio p they are 
all at rest (Quiescence (Q)) or all oscillating (Global Os- 
cillations (GO)). The region GO contains a thin stripe- 
shaped subregion, where all the units are oscillating, but 
with frequency ratios different from 1:1. 

The structure of the parameter space in Fig.^a) is re- 
produced for other oscillator types like the one proposed 
by Eguia et al. 15] and a phase model close to a SNIC 
(so-called 'active rotator' jig). The equations for the 
globally-coupled phase model are 

K N 

9 j = l- b : sin 0j + -J2 sin (^ ~ il ( 3 ) 
1=1 

where the parameter bj characterizes the non-uniformity 
of the phase rotation. Figure^b) shows that the simple 
phase model © qualitatively behaves as the ensemble 
of ML oscillators [Fig. Ufa)]. The only significant differ- 
ence is the shape of the region inside GO with different 
rotation numbers in each subpopulation: wedge-shaped 
for the phase model and stripe-shaped for "complete" 
models. For phase models (other implementations were 
tested), leaving the Q region implies generically a tran- 
sition to GO with 1 : 1 frequency ratio, or to PO (1 : 0). 
For coupled oscillators, the transition from the quiescent 
state to oscillations with other frequency ratios (1 : r) is 
of codimension one. From our simulations with the ML 
equations we have indications that in this transition the 
structure of phase space is a Cherry flow (see e.g. |17j), 
although this seems difficult to prove. 

In order to check the robustness of the diagrams in 
Fig. [21 we performed numerical simulations (N = 1000) 
considering a certain degree of heterogeneity in each sub- 
population (see red dotted lines in Fig. This has 
been carried out distributing the systems' parameters 
uniformly around 4>a,e (ML) and b A , e (phase model) 
with a finite width 7. Heterogeneity shrinks the region 
Q, since there is need of more coupling strength to com- 
pensate the major diversity in the population. As a last 
remark, we just notice that the dynamics in the PO and 
GO regions become more complex (some regimes cannot 
be labeled by just one frequency ratio). 

The rest of this Letter is devoted to the analysis of 
Fig. |21 in the absence of heterogeneity (i.e. 7 = 0). We 
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will mainly stress those results whose validity is general, 
i.e. irrespective of the particular dynamical system con- 
sidered. One empirical fact that simplifies the analysis is 
that in the whole K — p plane all the identical elements 
are at the same state: x je 5 4 = x^XjgSg = x^. Thus, 
it is possible to study the dynamics of the population in 
terms of two asymmetrically coupled elements: 

x A = F A (x A ) + Kp(x E - x A ), (4a) 
± E = F E {x E )+K(l-p)(x A -x E ), (4b) 

with x = (V, w) for the ML cells, and x = 9 for the phase 
model. This simplification allows to study the N — > oo 
limit and to adjust p continuously. In fact, reduction Q 
was used to efficiently compute Fig. [21 And in particu- 
lar for the phase model J3J), the bifurcation lines limit- 
ing Q can be easily calculated analytically for b A = 0: 
p c = l/K and p c — l/b El with a degenerate point at 
(K,p) = (b E , b^ 1 ). Other lines must be computed nu- 
merically EOj . 

An important feature of the phase diagrams in Fig. |3 
is the asymptotic value of the bifurcation line limiting Q: 
= p c (K — > oo). Remarkably, can be analytically 
estimated in a simple way. The calculation is based on 
the fact that close to p c the dynamics evolves for long 
time inside a (slow) region where the SN bifurcation takes 
place. Therefore, resorting to the normal form of a SN 
bifurcation, we obtain (after rcscaling space and time): 

z A = ae A - z\ + Kp(z E - z A ), (5a) 
z E = ae E - z% + K(l - p){z A - z E ), (5b) 

with e A < and e E > according to the oscillatory (no 
fixed point) and the excitable (one stable and one saddle 
fixed point) regimes at both sides of the SN bifurcation 
(at e = 0). In the limit K — > oo, the quadratic terms can 
be neglected at the bifurcation point. From the nullclines 
of Eq. (J3J) [z E = z A — (ae A — z\) / Kp and z A — z E — (ae E — 
z E )/K(l — p)) we obtain 



This is one of the main results of the present work. Equa- 
tion © gives the critical proportion of excitable and os- 
cillatory elements in order for the whole population to be- 
come self-oscillatory. Only three parameters are needed: 
the value of the control parameter at the SN bifurcation, 
and the distances, e A and e E , of the oscillatory and ex- 
citable elements to the bifurcation point. Specifically, for 
the ML model we have e A ^ E ~ {<p* — <p A . E ), that accord- 
ing to © yields p£° ~ 0.60, in good agreement with the 
numerical results (see also Fig. [3] below) 22]. Remark- 
ably we find that Eq. JfjJ) holds also for an ensemble of 
dynamical systems at both sides of a -ffop/bifurcation (we 
skip here the proof, but cf. Eq. (4) in [j] for a particular 
case). 

Our next results concern the behavior of the mean en- 
semble's frequency, which is a natural measure for the 
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FIG. 3: (a) Normalized ensemble's frequency 0,(p) for dif- 
ferent values of the coupling strength K in a population of 
oscillatory and excitable ML units (TV = 500) . The dashed 
line for large coupling (K — 5) shows the convergence of p c 
to the value predicted by Eq. ©. Panels (b-e) show the in- 
dividual frequencies for: K = 0.05 (PO)(top), K c ~ 0.144, 
K = 0.2, and K = 0.3 (transition GO-Q) (bottom). 

global oscillatory activity. The ensemble's average of the 
individual frequencies is 

1 N 

3=1 

where ttj represents the mean frequency of the j-th el- 
ement. Obviously, under the assumption in Eq. |0J, 
Q = (1 -p)Q A +pVt E (with Vl A = Q-jeSAi = &jeS E )- 
Fig. |3 shows the average frequency Q — normalized by 
the frequency at p — 0: fl(p) = O(p)/O(0) — , and the 
individual frequencies Qj(p) (right panels) for different 
values of K. Notice that for small coupling the transition 
to Quiescence (f2 = 0) occurs trivially at p c = 1. Such 
behavior changes above K c ~ 0.144, where the transi- 
tion begins to take place at p c < 1. Also, note that 
for intermediate values of K the curves present a step- 
like profile, corresponding to the stripe-shaped region in 
Fig. |2a). This can be seen in Figs.|2Ic,d) where, in the 
corresponding range of p, the two observed frequencies 
in the ensemble arc not 1 : 1 (neither 1 : 0) related. 

The quantity f2 presents interesting universal proper- 
ties around K c . Assuming that Eq. exhibits a SN 
bifurcation at some p sn (K) we have : p c — 1 < p sn (K) 
(for K < K c ), p c = 1 = p sn (K) (at K = K c ) and 
Pc = Psn (K) < 1 (for K > K c ) . And noting that 

(i) Only the active (self-oscillating) elements con- 
tribute to ft: tt = (1 -p)fl A . [Tl E — 0] . 

(ii) In a SNIC, the frequency of the oscillating elements 
scales as a square root: fl A ~ h(p sn — p) 1 / 2 . 

we find that £1 grows from zero as a power of the distance 
to the critical p c : 

n = (i- p)h{ Psn - p) 1 ' 2 cx (p c - P f, (8) 
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FIG. 4: Fitting to Eq. JTOJ for an ensemble of TV = 1000 ML 
elements, K c = 0.14386(0). Dashed lines stem from Eq. itTUt 
with fitting parameters g = 5.01 and h = 0.66. 



with three different values of the exponent (3: 1 (for K < 
K c , Psn > V c = 1), 3/2 (at K = K c , p sn = p c = 1), 1/2 
(for K > K c , p sn =p c <l). 

In a neighborhood of (K = K c ,p = 1), we may obtain 
a universal scaling function by assuming that the shift of 
p sn is approximately linear on K: p sn ~ 1 + g{K c — K). 
Recalling that p c — 1 for K < K c , we obtain from Eq. (JSJ) 
the expressions 



n = h(p c ~p)[p c -p + g(K c -K)}^ 2 , (K<K C ) 

,1/2 (V ^ TS \ (J) 



n = h(i - p)( Pc - P y 



(K > K c ) 



These distinct scalings may be condensed into the single 
scaling function (by approximating 1 ~ p c — g(K c — K) 
in the second equation): 



n = h{p c ~pf' 2 <S> (g 



K-K c 
Pc-P 



(10) 



with ®(x) [= v 7 ] -7 ~x (x < 0)], [= 1 +x(x > 0)]. We 
obtain then two scaling regions with common fitting pa- 
rameters g and h, see Fig. 0J 

It is important to note that the scaling in Fig. 0] co- 
incides with the one reported in Q for the amplitude of 
oscillations in a mixed population of oscillators close to 
a Hopf bifurcation. Such coincidence lies in the same 
asymptotic dependence (square-root law) for the cycle's 
amplitude in the case of a Hopf bifurcation and for the 
cycle's frequency in the case of a SNIC. 

In conclusion, we have demonstrated that for systems 
close to a SNIC, a transition to global quiescence oc- 
curs as the ratio p of excitable elements in the ensemble 
exceeds a certain value p c . Such transition is a gener- 
alization of the Aging Transition reported in £j and is 
characterized by a universal scaling function relating the 
mean frequency Q with p and the coupling strength K. 
Additionally, we derive an analytical estimation for p c in 
the large K limit which holds for both Hopf- and SNIC- 
mediated Aging Transitions. 

Our results might be of importance in several situa- 
tions, since excitability is a typical feature of many phys- 
ical, chemical and biological systems. Hence, our work 
is a first step in modeling the competition between ex- 
citable and oscillatory dynamics, with possible extensions 
to extended media and complex networks. Finally, our 
findings could also be relevant in populations of excitable 
systems when some of the elements turn self-oscillatory 
due to the presence of noise (coherence resonance ll8j). 
a situation common in neuroscience. 
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G. Kalosakas for a critical reading of the manuscript. 
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